Overview

The role of temporal dispersal patterns in building diverse tallgrass prairie plant communities

Background:

Establishing diverse plant communities is critical since diversity is linked to ecosystem health. However, recreating tallgrass prairies with high plant diversity has been challenging, especially for early-season species. Variation in species arrival may influence plant community composition and diversity. Since prairies exhibit temporal patterns of seed dispersal, planting all species simultaneously could forgo phenological differences that promote species coexistence. Therefore, we investigated whether manipulating plant species’ arrival according to natural dispersal phenology influences reconstruction outcomes.

In 2021, we manipulated the arrival of 36 species via seed additions of (i) species in the order of peak dispersal timing, (ii) summer dispersing species (first peak in dispersal activity before September 1st) followed by 18 fall dispersing species, (iii) 18 fall dispersing species followed by 18 summer dispersing species, (iv) all species simultaneously. Additionally, we had a negative control that had no seed additions. One year later, we found that differences in seeding treatment influenced the diversity and composition of reconstructed communities. Species arriving later had less cover than when seeded with priority, particularly for summer dispersing species. Overall, our study provides evidence of priority effects in reconstructed grasslands and suggests that the timing of seed additions post-disturbance influences restoration outcomes.

Author: Katherine Carter Wynne ()

Co-authors: Lauren L. Sullivan

Created: 06 December 2022


Files:

  1. PE_Cover_Inner_Fall_2021_Cleaned.xlsx - vegetative cover taken at peak biomass in 2021 (August 27th - August 30th); only half of the lumped seedings occured.

  2. PE_Cover_Inner_Summer_2022_Cleaned.xlsx - vegetative cover measured in Summer 2022 (June 27th - June 28th)

  3. PE_Cover_Inner_Fall_2022_Cleaned.xlsx - vegetative cover taken at peak biomass in 2022 (August 25th - August 26th)

  4. PE_Cover_Inner_Summer_2023_Cleaned.xlsx - vegetative cover measured in Summer 2023 (June 21st - June 23rd)

  5. PE_Cover_Inner_Fall_2023_Cleaned.xlsx - vegetative cover taken at peak biomass in 2023 (August 10th - August 11th)

  6. Bradford_10_Year_Weather.xlsx - weather (temperature and precipitation) data from the last 10 years at Bradford Research Farm

  7. Priority_Effects_Fall_2022_Light.xlsx - light data from August 2022; measured at the same time as peak biomass cover

  8. Priority_Effects_Fall_2023_Light.xlsx - light data from August 2023; measured at the same time as peak biomass cover

  9. PE_Bare_Ground_Cleaned_2021.xlsx - Litter and bare ground cover from August 2021

  10. Species_List.xlsx - Species list


R Version: R 4.2.2

RStudio Version: 2023.06.1+524

Package Version: Found in the .readme document in GitHub Repository

Last updated: 29 September 2023


Overall Setup

Materials and Methods

Study Site

Our experiment utilized a former agricultural field (~1 acre) at the University of Missouri’s Bradford Research Center in Columbia, MO (38.893604, -92.201154, Boone County, MO). Typical for the central U.S., the 10-year (2011 – 2021) mean annual precipitation and average air temperature at Bradford Research Center were 927.74 \(\pm\) 143.01 mm and 12.40 \(\pm\) 10.45 \(\circ\)C, respectively (Commercial Agriculture Automated Weather Station Network, d.n.). Prior to our experiment, the field we grew herbicide-resistant soybeans for at least three years, reflecting conditions similar to most prairie reconstructions before seeding (Rowe, 2010; Newbold et al., 2019). In March 2021, we tilled our study site and hand-removed rhizome clumps to create a smooth surface before our first seeding.

Experimental Design

We used seed additions to manipulate the arrival of 36 native tallgrass prairie plant species. Since prairies experience two peaks in dispersal activity (Rabinowitz & Rapp 1980), we classified species into two dispersal guilds, summer-dispersing species (first peak in dispersal activity before September 1st) and fall-dispersing species. Dispersal guild is more informative than flowering guild because dispersal does not always follow flowering (e.g., Penstemon digitalis). We based our classifications on expert opinion from the Missouri Department of Conservation, The Nature Conservancy, and our previous work on seed rain patterns in Missouri tallgrass prairies (Wynne et al. in prep). Species used in our study consisted of 29 native species captured in our previous study on seed rain and seven summer-dispersing species that restoration managers cited as having minimal success in prairie reconstructions (e.g., Viola pedatifida) (Table S1) (Barak et al., 2022). When possible, we obtained seeds from local ecotype commercial sellers. However, several summer-dispersing species were not grown commercially in Missouri and were sourced elsewhere. We stored the seeds in a refrigerator (2.78 \(\circ\)C) until seeding.

We used a randomized factorial block design to test the effects of seeding timing and order (Fig. 1). Each block (n = 6) contained five plots (2 x 2 m) randomly assigned a seed addition treatment following one of the four arrival treatments:

  1. the addition of species in order of first peak in dispersal activity (NAT)

    - May 24th 2021 (VIOPED, PACPLA)
    - June 4th 2021 (SISCAM, SPHOBT, CORLAN)
    - June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC)
    - July 11th 2021 (CORPAL, DODMEA)
    - July 25th 2021 (KOEMAC, AMOCAN)
    - August 2nd 2021 (LINSUL, MELVIR)
    - August 21st 2021 (DALCAN, DALPUR, ACHMIL)
    - September 17th 2021 (CROSAG)
    - October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR)
    - October 18th 2021 (BIDARI, HELMOL, LESCAP)
    - November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC)
    - November 19th 2021 (CORTRI, PYCTEN, SOLRIG)
  2. the lumped addition of 18 summer-dispersing species on March 22nd, 2021 followed by a lumped addition of 18 fall-dispersing species approximately five months later on September 5th, 2021 (LE)

  3. the lumped addition of 18 fall-dispersing species on March 22nd, 2021 followed by a lumped addition of 18 summer-dispersing species approximately five months later on September 5th, 2021 (LL)

  4. the simultaneous addition of the entire species pool on March 22nd, 2021 typical of most prairie reconstructions (SIM).

Additionally, we had an unseeded negative control (NON) to determine whether active seeding treatments improved community diversity compared to passive regeneration from ambient seed sources (e.g., seed rain and soil seed bank). We also seeded white clover (Trifolium repens) between the experimental plots to prevent erosion. In fall 2022, we started annually removing invading white clover and woody species from plots. All other unseeded species were not weeded to better reflect realistic assembly and reconstruction processes.

We started seeding at the end of the dormant season (March 22nd, 2021) and continued until November 19th, 2021. We hand-seeded species into experimental plots at a density of 50 seeds \(m^{-2}\) using sand as a broadcasting agent. For treatments requiring multiple seedings, we incorporated M-Binder tackifier (Ecology Controls, Carpinteria, CA), a natural adhesive often used in hydroseeding, into the seeding mixes to increase soil-seed contact without disturbing the existing vegetation. After seeding, we lightly watered plots to activate the tackifier. All plots received equal amounts of sand, tackifier, and water during every seeding.

Data Collection

Starting in 2022, we conducted floristic surveys to assess plant community diversity and composition in experimental plots for two years. To determine species abundance, we measured the percent aerial cover of all vascular plant species rooted within a 1 \(m^{2}\) subplot in the center of the experimental plot. Because many seeded summer-dispersing species senesce before peak biomass (August – September), we measured vegetative cover twice yearly, once in the early summer (late June) and again at peak biomass (late August). We only used the highest percent cover value for species present in both surveys. We identified species according to Yatskievych (1999, 2006, 2013). Additionally, we measured environmental variables including percent bare ground cover and light levels in plots. For each plot, we measured light levels using a (insert name of light bar) two times above the vegetation, in the midstory, and ground-level at solar noon \(\pm\) 2 hours.

Data Setup

Each tab below contains the code used to import, clean, and manipulate data to later use for further analyses.

Load libraries

## Libraries


### ----- Data Cleaning and Management -----

# Import excel files

library(readxl)


# Export dataframes to excel files

library(writexl)


# Data cleaning and visualization

library(tidyverse)



### ----- Data Analysis -----

# Community ecology functions (NMDS, diversity indices, etc.)

library(vegan)
library(BiodiversityR)

# Mixed models

library(lme4)
library(lmerTest)
library(MuMIn)

# Checks model assumptions

library(rstatix)

### ----- Data visualization -----


#Functions that assist in making high quality figures

library(cowplot)
library(ggpubr)
library(ggrepel)
library(flextable)

Import Datasets

## Import Datasets

## inner (1 m^2) cover for Summer
#### comprehensive vegetative cover (%) for all species in the inner (1 m^2) permanent sampling area. 

inner_cover_summer_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Summer_2022_Cleaned.xlsx")

inner_cover_summer_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Summer_2023_Cleaned.xlsx")

## inner (1 m^2) cover for Fall
#### comprehensive vegetative cover (%) for all species in the inner (1 m^2) permanent sampling area. 


inner_cover_fall_2021 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2021_Cleaned.xlsx")

inner_cover_fall_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2022_Cleaned.xlsx")

inner_cover_fall_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Cover_Inner_Fall_2023_Cleaned.xlsx")


## Bradford weather
weather <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Bradford_10_Year_Weather.xlsx")


## Light

light_fall_2022 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Priority_Effects_Fall_2022_Light.xlsx")

light_fall_2023 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Environment/Priority_Effects_Fall_2023_Light.xlsx")


## Bare and litter - 2021

bare_cover_2021 <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Cover/PE_Bare_Ground_Cleaned_2021.xlsx")


## Species list

species_list_PE <- read_excel("~/Priority-Effects/Priority Effects - Github/Data/Species_List.xlsx")

Dataset cleaning

2021

## 2021 Inner Cover - Dataset cleaning


### Remove unnecessary columns

#### Get rid of the notes, unknown, and senensced columns

### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]


#### change log
 # - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct.  (Sep 6, 2023)
 # - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)


### Get rid of EMP

### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")

2022

## 2022 Inner Cover - Dataset cleaning

### Remove unnecessary columns

### Get rid of the notes, unknown, and sentenced columns

#### Summer
inner_cover_summer_2022 <- inner_cover_summer_2022[,-8]

#### change log
 # - Fixed IDs: most GALAPA -> GALPED and HEURIC -> GEUCAN (Sep 6, 2023)


#### Fall
inner_cover_fall_2022 <- inner_cover_fall_2022[,-c(8,9,10)]

#### change log
 # - Fixed IDs: VERSPP -> VERARV since this is the species observed in the summer (Sep 6, 2023)



### Get rid of EMP 

#### Summer
inner_cover_summer_2022 <- subset(inner_cover_summer_2022, Treatment != "EMP")

#### Fall
inner_cover_fall_2022 <- subset(inner_cover_fall_2022, Treatment != "EMP")

2023

## 2023 Inner Cover - Dataset cleaning

### Remove unnecessary columns

### Get rid of the notes, unknown, and sentenced columns

#### Summer
inner_cover_summer_2023 <- inner_cover_summer_2023[,-c(6,7)]

#### change log
 # - Fixed IDs: was able to ID ALLSPP as ALLVIN (Sep 6, 2023)


#### Fall
inner_cover_fall_2023 <- inner_cover_fall_2023[,-c(6,7)]


### Get rid of EMP 

#### Summer
inner_cover_summer_2023 <- subset(inner_cover_summer_2023, Treatment != "EMP")

#### Fall
inner_cover_fall_2023 <- subset(inner_cover_fall_2023, Treatment != "EMP")

Joining Years

#Checking to make sure all the SPP6 levels are correct (no misspellings or weirdness)

### Fall 2021
sort(unique(inner_cover_fall_2021$SPP6))
##  [1] "ACAVIR"     "ACHMIL"     "AGRHYE"     "AMATUB"     "AMBART"    
##  [6] "AMBTRI"     "BARVUL"     "BIDARI"     "BROJAP"     "CARFRA"    
## [11] "CHAFAS"     "CONCAN"     "CORLAN"     "CORTIN"     "CORTRI"    
## [16] "CROSAG"     "CYNDAC"     "CYPODO"     "DALCAN"     "DIGISC"    
## [21] "DIGSAN"     "ECHSPP"     "ERACIL"     "ERIANN"     "ERYYUC"    
## [26] "GERCAR"     "GLYMAX"     "HELAUT"     "HELMOL"     "HORPUS"    
## [31] "IPOHED"     "IPOLAC"     "IVAANN"     "JUNINT"     "LEPSPP"    
## [36] "LESCAP"     "LINSUL"     "LONMAC"     "LOTCOR"     "MELOFF"    
## [41] "MOLVER"     "MONFIS"     "MYOVER"     "OXADIL"     "PACKERA?"  
## [46] "PANCAP"     "PANDIC"     "PENDIG"     "PERHYD/PUN" "PERPEN"    
## [51] "PHAARU"     "PHYSPP"     "PLALAN"     "PLAMAJ"     "PLAVIR"    
## [56] "POTNOR"     "RANARV"     "RATPIN"     "RUDHIR"     "RUMCRI"    
## [61] "SCHSCO"     "SETFAB"     "SETPUM"     "SIDSPI"     "SOLCAN"    
## [66] "SOLCAR"     "SORNUT"     "SYMPIL"     "TAROFF"     "TRIREP"    
## [71] "VERURT"     "VIOPED"     "XANSTR"
### Summer 2022
sort(unique(inner_cover_summer_2022$SPP6))
##  [1] "ACHMIL" "AGRHYE" "BARVUL" "BROJAP" "CARBRE" "CARFRA" "CERSPP" "CORLAN"
##  [9] "CORPAL" "ERIANN" "GALAPA" "GALPED" "GERCAR" "GEUCAN" "HORPUS" "LOTCOR"
## [17] "MELOFF" "MYOVER" "PENDIG" "PHAARU" "PLALAN" "PLAMAJ" "PLAVIR" "POTNOR"
## [25] "RANABO" "RUDHIR" "RUMCRI" "SILANT" "SPHOBT" "TAROFF" "TRIPRA" "TRIREP"
## [33] "VERARV" "VIOPED"
### Fall 2022
sort(unique(inner_cover_fall_2022$SPP6))
##  [1] "ACAVIR"     "ACHMIL"     "AGRHYE"     "AMBART"     "AMBTRI"    
##  [6] "BARE"       "BARVUL"     "BIDARI"     "BROJAP"     "CARFRA"    
## [11] "CARSPP"     "CHAFAS"     "CONCAN"     "CORLAN"     "CORPAL"    
## [16] "CORTRI"     "CYNDAC"     "DESSPP"     "DIGISC"     "DIGSAN"    
## [21] "ECHSPP"     "ERASPE"     "ERIANN"     "ERYYUC"     "FESARU"    
## [26] "HELAUT"     "HELMOL"     "HORPUS"     "IPOLAC"     "IVAANN"    
## [31] "JUNINT"     "LACSER"     "LESCAP"     "LIAPYC"     "LINSUL"    
## [36] "LITTER"     "LONMAC"     "LOTCOR"     "MONFIS"     "MORALB"    
## [41] "MYOVER"     "OXADIL"     "PANCAP"     "PANDIC"     "PENDIG"    
## [46] "PERHYD/PUN" "PHAARU"     "PLALAN"     "PLAMAJ"     "POAPRA"    
## [51] "POTNOR"     "RANABO"     "RATPIN"     "RUDHIR"     "RUMCRI"    
## [56] "SETFAB"     "SETPUM"     "SIDSPI"     "SOLCAN"     "SOLCAR"    
## [61] "SORNUT"     "SPOHET"     "TAROFF"     "TRIPRA"     "TRIREP"    
## [66] "VERARV"     "VERURT"     "VIOPED"
### Summer 2023
sort(unique(inner_cover_summer_2023$SPP6))
##  [1] "ACHMIL" "AGRHYE" "ALLVIN" "BARVUL" "BROJAP" "CARBRE" "CARFRA" "CARSPP"
##  [9] "CERSPP" "CORLAN" "ERIANN" "FESARU" "GALAPA" "GALPED" "GERCAR" "GEUCAN"
## [17] "HORPUS" "JUNINT" "KOEMAC" "LEPVIR" "LONMAC" "LOTCOR" "MELSPP" "MYOVER"
## [25] "OXADIL" "PENDIG" "PHAARU" "PLALAN" "PLAMAJ" "PLAVIR" "POAPRA" "PYCTEN"
## [33] "RANABO" "RUDHIR" "RUMCRI" "SILANT" "SPHOBT" "TAROFF" "TRIPRA" "TRIREP"
## [41] "VERARV" "VIOPED"
### Fall 2023
sort(unique(inner_cover_fall_2023$SPP6))
##  [1] "ACHMIL" "AGRHYE" "AMBART" "AMBTRI" "BARE"   "BARVUL" "BIDARI" "CARFRA"
##  [9] "CARSPP" "CERSPP" "CHAFAS" "CONCAN" "CORLAN" "CORTRI" "CYNDAC" "DESSPP"
## [17] "ECHSPP" "ERASPE" "ERIANN" "ERYYUC" "FESARU" "GEUCAN" "HELAUT" "HELMOL"
## [25] "IPOHED" "IPOLAC" "IVAANN" "KOEMAC" "LESCAP" "LIAPYC" "LITTER" "LONMAC"
## [33] "LOTCOR" "MELALB" "MELSPP" "MONFIS" "MYOVER" "OXADIL" "PANDIC" "PENDIG"
## [41] "PHAARU" "PLALAN" "PLAMAJ" "POAPRA" "PYCTEN" "RANABO" "RATPIN" "RUDHIR"
## [49] "RUMCRI" "SETFAB" "SETPUM" "SIDSPI" "SOLCAN" "SOLCAR" "SORNUT" "SPOHET"
## [57] "SYMPIL" "TAROFF" "TRIPRA" "TRIREP" "VERURT" "VIOPED"
### ^ above looks good so far
# Make a unique identifier for year and season

### 2021

#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))


### 2022

# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)

### 2023

#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))

#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)

inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)
# Join summer and fall cover for each year 

full_2021_data <- inner_cover_fall_2021

full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)

full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)


# Join all years together

### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)

### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)
# lump certain taxa together 

### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)

## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP

for(i in 1:nrow(full_year_data)) {
  if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
  if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
  if(full_year_data[i,3] == "LEPSPP"){full_year_data [i,3] <- "LEPVIR"}
  if(full_year_data[i,3] == "ECHSPP"){full_year_data [i,3] <- "ECHCRU"}
  
  # Make all cover data that was less than 1 = to 1 (to make the data discrete)
  if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
  
}
# Function below takes that highest cover value for a species seen twice (e.g., once in the summer and again in fall)

inner_cover_max <- full_year_data  %>%
      group_by(Block, Treatment, Year, SPP6) %>%
      summarise(max_cover=max(Percent_Cover))

##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice


### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")


### Remove Litter
inner_cover_max_only  <- subset(inner_cover_max_only, SPP6 != "LITTER")
#Checking again to make sure all the SPP6 levels are correct (no misspellings or weirdness)

sort(unique(inner_cover_max_only$SPP6))
##  [1] "ACAVIR"     "ACHMIL"     "AGRHYE"     "ALLVIN"     "AMATUB"    
##  [6] "AMBART"     "AMBTRI"     "BARVUL"     "BIDARI"     "BROJAP"    
## [11] "CARSPP"     "CERSPP"     "CHAFAS"     "CONCAN"     "CORLAN"    
## [16] "CORPAL"     "CORTIN"     "CORTRI"     "CROSAG"     "CYNDAC"    
## [21] "CYPODO"     "DALCAN"     "DESSPP"     "DIGISC"     "DIGSAN"    
## [26] "ECHCRU"     "ERACIL"     "ERASPE"     "ERIANN"     "ERYYUC"    
## [31] "FESARU"     "GALAPA"     "GALPED"     "GERCAR"     "GEUCAN"    
## [36] "GLYMAX"     "HELAUT"     "HELMOL"     "HORPUS"     "IPOHED"    
## [41] "IPOLAC"     "IVAANN"     "JUNINT"     "KOEMAC"     "LACSER"    
## [46] "LEPVIR"     "LESCAP"     "LIAPYC"     "LINSUL"     "LONMAC"    
## [51] "LOTCOR"     "MELSPP"     "MOLVER"     "MONFIS"     "MORALB"    
## [56] "MYOVER"     "OXADIL"     "PACKERA?"   "PANCAP"     "PANDIC"    
## [61] "PENDIG"     "PERHYD/PUN" "PERPEN"     "PHAARU"     "PHYSPP"    
## [66] "PLALAN"     "PLAMAJ"     "PLAVIR"     "POAPRA"     "POTNOR"    
## [71] "PYCTEN"     "RANABO"     "RANARV"     "RATPIN"     "RUDHIR"    
## [76] "RUMCRI"     "SCHSCO"     "SETFAB"     "SETPUM"     "SIDSPI"    
## [81] "SILANT"     "SOLCAN"     "SOLCAR"     "SORNUT"     "SPHOBT"    
## [86] "SPOHET"     "SYMPIL"     "TAROFF"     "TRIPRA"     "TRIREP"    
## [91] "VERARV"     "VERURT"     "VIOPED"     "XANSTR"
SPP6_list <- sort(unique(inner_cover_max_only$SPP6))
SPP6_list <- as.data.frame(SPP6_list)

### Saw at least 94 different species across study
#### More since I had to lump MELOFF, MELALB, CARFRA, and CARBRE


### Function below makes the species list 

# write_xlsx(SPP6_list, "Species_List2.xlsx")
### Add species information to the dataset

#### N/A is N = native, A = non-native, G = Genus

inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)

#Remove unnecessary columns

inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]

Bare/Litter

### Bare dataset

inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")

### Litter dataset

inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")


### Manipulate 2021 dataset to match 2022 and 2023

#### Bare Ground

bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))

for(i in 1:nrow(bare_2021)) {
  if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}

#### Litter 

litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))


for(i in 1:nrow(litter_2021)) {
  if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}


### Full_join datasets

inner_bare_cover <- full_join(inner_bare_cover, bare_2021)

inner_litter_cover <- full_join(inner_litter_cover, litter_2021)

Question 1: Does timing of species arrival influence local diversity?

Setup and Data Management

Data Analysis

We fit mixed-effects linear models with block as a random effect to determine whether year, timing of species arrival via seeding treatments, and the interaction between year and seeding treatments influenced local diversity (i.e., species richness, mean conservatism value, and Shannon diversity index) in reconstructed tallgrass prairie plant communities. In cases where random effect variance was estimated as near zero, we dropped the random effect and fit a simpler linear model instead. For each diversity model, we conducted a type III analysis of variance (ANOVA) followed by a posthoc Tukey test to determine significant pairwise differences.

Below is the code used to calculate the species richness, mean conservatism value, and Shannon diversity index for each plot, treatment, and year combination.

### Making a summary table of diversity indices between treatments

#Make a unique identifier for every block and treatment
inner_cover_max_reduced$UniqueBlock <- paste(inner_cover_max_reduced$Treatment, inner_cover_max_reduced$Block, sep="_")


### Making a community matrix

inner_cover_max_reduced_lumped <- inner_cover_max_reduced %>% 
  filter(Year != "2021") %>% 
  group_by(Block, Treatment, Year, SPP6) %>% 
  summarise(max_cover = sum(max_cover))


#Make a wide formatted dataset

inner_wide <- inner_cover_max_reduced_lumped  %>% 
  spread(key="SPP6", value="max_cover") 

#Replace NA values with 0

inner_wide[is.na(inner_wide)] <- 0 

#Make a separate dataframe for the labels
inner_wide.labs <- inner_wide[, c(1,2,3)]

#Turn our dataset into a matrix 
inner_wide_mat <- inner_wide[,-c(1,2,3)]


#Double check your work using the View() function. 

#View(inner_wide_mat)
## Calculate diversity metrics


# --- Species richness --- 

Species_richness <- specnumber(inner_wide_mat)


#  --- Shannon's Diversity --- 

#Calculate Shannon's Diversity
Shannon <- diversity(inner_wide_mat, index="shannon")


# --- Mean C --- 


# Make a dataset that has each treatment, the species found in that treatment, and their associated C values

inner_cover_max_reduced$C_Value <- as.numeric(inner_cover_max_reduced$C_Value)

C_hist <- inner_cover_max_reduced  %>% 
                  filter( C_Value >= 0) %>% 
                  group_by(Block, Treatment, SPP6, Year, C_Value) %>% 
            summarize(tot.cover = sum(max_cover)) %>% 
            select(-c(tot.cover))


# Calculate mean C for each Treatment in a Block
### Remember you have to filter out anything that isn't native!

Only_C <- inner_cover_max_reduced  %>% 
                  filter(Year != "2021") %>% 
                  filter( C_Value >= 0) %>% 
                  filter( C_Value != "NA") 

Only_C$C_Value <- as.numeric(Only_C$C_Value)

Mean_C <- Only_C %>% 
                  group_by(Block, Treatment, Year) %>% 
                  summarize(Mean_C = mean(C_Value))

Diversity Indices - Summary Statistics

Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.

Treatment

Year

C

± SD

Richness

± SD

SDI

± SD

LE

2022

1.68

0.58

20.00

3.46

2.12

0.20

LE

2023

1.91

0.47

18.17

1.83

2.18

0.16

LL

2022

1.99

0.41

21.00

2.19

2.25

0.20

LL

2023

2.62

0.20

17.83

2.40

2.16

0.13

NAT

2022

1.43

0.84

17.67

2.34

2.05

0.12

NAT

2023

1.45

0.73

17.00

3.69

2.01

0.28

NON

2022

0.48

0.34

15.00

2.61

1.86

0.17

NON

2023

0.42

0.41

14.50

2.74

1.75

0.20

SIM

2022

2.34

0.44

24.50

3.89

2.23

0.32

SIM

2023

2.84

0.32

22.17

5.08

2.34

0.36

Mean C

BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 16, family="Arial"),
axis.text = element_text(size = 13, colour = "gray25"),
axis.title = element_text(size = 16, colour = "gray25"),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16),
legend.key = element_blank())

Mean C - Model assumptions

Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. Simultaneous seeding had the most “even” distribution of C values and almost every value was represented in this treatment. Across years, the largest changes in the distribution of C values were in the treatments that seeded the fall dispersing species in September. These changes were attributed to the replacement of species with low c values with more conservative species. Plots that were never seeded mainly consisted of species with low conservatism value.

Generally, mean C values were normally distributed with few extreme outliers, suggesting we can use a linear model to predict mean c as a function of treatment and year.

Mean C - Model fitting

I used a linear model to predict mean c as a function of seeding treatment, year, and an interaction between treatment and year. A mixed model using block as a random effect produced a singular fit because the variance explained by the random effect block was estimated as close to zero and did not further inform the data. Therefore, I chose to use the simpler model.

Model fit was high, with an adjusted \(R^{2}\) value of 0.687. Also, residual plots did not reveal any concerning patterns.

## 
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.07124 -0.32383 -0.07948  0.28266  1.12107 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.6602     0.1608  10.327 2.16e-14 ***
## TreatmentLL    0.5149     0.2075   2.481 0.016251 *  
## TreatmentNAT  -0.3549     0.2075  -1.710 0.093003 .  
## TreatmentNON  -1.3445     0.2075  -6.479 2.89e-08 ***
## TreatmentSIM   0.7983     0.2075   3.847 0.000318 ***
## Year2023       0.2660     0.1313   2.026 0.047691 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5084 on 54 degrees of freedom
## Multiple R-squared:  0.7134, Adjusted R-squared:  0.6869 
## F-statistic: 26.88 on 5 and 54 DF,  p-value: 1.544e-13
# Residual plot looks decent
plot(Mean_C.mod.lm.additive)

Analysis of variance (type III) revealed that seeding treatment was the only significant factor influence mean c value.

## Anova Table (Type III tests)
## 
## Response: Mean_C
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 27.561  1 106.6537 2.160e-14 ***
## Treatment   33.674  4  32.5772 8.088e-14 ***
## Year         1.061  1   4.1057   0.04769 *  
## Residuals   13.955 54                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc analysis indicated that all seeding treatments increased mean c value compared to unseeded plots. Seeding all species in one addition produced more conservative communities than treatments that seeded summer dispersing species first followed by fall dispersing species. Even though seeding only the fall dispersing species first produced higher quality communities than those that seeded species with natural timing, there was no difference in mean c between the seeding treatments that had two lumped additions.

Mean_C.mod.lm.est <- lm(Mean_C~Treatment, data =Summary_table)


est.lm.treatment_treatment <- emmeans::emmeans(Mean_C.mod.lm.est, ~ Treatment, type = "response")

pairs(est.lm.treatment_treatment, adjust = "tukey")
##  contrast  estimate    SE df t.ratio p.value
##  LE - LL     -0.515 0.213 55  -2.414  0.1270
##  LE - NAT     0.355 0.213 55   1.664  0.4644
##  LE - NON     1.345 0.213 55   6.303  <.0001
##  LE - SIM    -0.798 0.213 55  -3.742  0.0039
##  LL - NAT     0.870 0.213 55   4.077  0.0013
##  LL - NON     1.859 0.213 55   8.717  <.0001
##  LL - SIM    -0.283 0.213 55  -1.329  0.6748
##  NAT - NON    0.990 0.213 55   4.639  0.0002
##  NAT - SIM   -1.153 0.213 55  -5.406  <.0001
##  NON - SIM   -2.143 0.213 55 -10.045  <.0001
## 
## P value adjustment: tukey method for comparing a family of 5 estimates
#Posthoc test
#est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.interaction, ~ Treatment*Year, type = "response")

#pairs(est.lm.treatment, adjust = "tukey")

Mean C - Figures

Overall, seeding fall dispersing species early either as part of one simultaneous or as two lumped additions produced the highest quality communities.

Species Richness

SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.

SR - Model assumptions

Overall, species richness follows a normal distribution despite being count data. There were a couple of outliers in the natural seeding treatment but this is likely alright.

SR - Model fitting

I chose to use a linear mixed-effects model with block as a random effect to predict species richness as a function of seeding treatment, year, and an interaction between treatment and year. Despite being count data, the distribution of species richness is remarkably balanced and not skewed. Model comparison using AIC revealed that the Normal model performed better than the Poisson model, which was underdispersed. A more complicated Quasipoisson model that accounted for underdispersion produced similar results to the Normal mixed model. Therefore, I chose to use the less complicated Normal model.

Overall, model fit was decent. For example, \(R^{2}_{GLMM(m)}\) (i.e., the variance explained by only the fixed effects) was 0.474 and \(R^{2}_{GLMM(c)}\) (i.e., variance explained by the entire model) was 0.598. Inclusion of the random effect block substantially improved model fit. Furthermore, there were no strange patterns or occurrences in the model residuals.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 231.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0095 -0.7513  0.1272  0.6940  1.5317 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 1.828    1.352   
##  Residual             8.513    2.918   
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   18.3333     1.0915 23.7757  16.796 1.09e-14 ***
## TreatmentLE    1.7500     1.1912 38.0000   1.469   0.1500    
## TreatmentLL    2.0833     1.1912 38.0000   1.749   0.0884 .  
## TreatmentSIM   6.0000     1.1912 38.0000   5.037 1.18e-05 ***
## Year2023      -2.0000     0.8423 38.0000  -2.375   0.0227 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLE TrtmLL TrtSIM
## TreatmentLE -0.546                     
## TreatmentLL -0.546  0.500              
## TreatmntSIM -0.546  0.500  0.500       
## Year2023    -0.386  0.000  0.000  0.000
##            R2m       R2c
## [1,] 0.3644925 0.4768066

A type III ANOVA revealed significant differences in species richness between years and seeding treatments. However, there was no significant interactive effect between year and seeding treatment on species richness.

Seeding treatments that added species either in one or two lumped additions had significantly greater species richness than unseeded plots. Multiple seedings according to natural dispersal activity did not increase species richness compared to unseeded plots. Adding the entire species pool simultaneously produced communities with the greatest species richness. Treatments that conducted multiple seedings resulted in communities possessing similar amounts of species richness regardless of timing and order of species arrival.

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Species_richness
##                Chisq Df Pr(>Chisq)    
## (Intercept) 282.1205  1  < 2.2e-16 ***
## Treatment    27.1051  3  5.596e-06 ***
## Year          5.6383  1    0.01757 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast  estimate   SE df t.ratio p.value
##  NAT - LE    -1.750 1.26 39  -1.389  0.5137
##  NAT - LL    -2.083 1.26 39  -1.653  0.3616
##  NAT - SIM   -6.000 1.26 39  -4.762  0.0002
##  LE - LL     -0.333 1.26 39  -0.265  0.9934
##  LE - SIM    -4.250 1.26 39  -3.373  0.0088
##  LL - SIM    -3.917 1.26 39  -3.108  0.0177
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

SR - Figures

Overall, seeding all species at the same time produced communities with the greatest species richness. For the most part, seeding increased species richness

Summer-Dispersing

### Filter out everything that isn't a seeded summer-dispersing species

inner_long <- inner_wide %>% 
  gather(key = "SPP6", value = "max_cover", ACAVIR:VIOPED) 

Sown_species_ID <- left_join(inner_long,  species_list_PE)
## Joining with `by = join_by(SPP6)`
SR_summer <- Sown_species_ID  %>% 
  group_by(Block, Treatment, Year) %>% 
  filter(Seeded == "Summer" & max_cover > 0 & Treatment != "NON") %>% 
  summarize(Summer_SR = n_distinct(SPP6))
## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.
SR_summer_zero <- Sown_species_ID  %>% 
  group_by(Block, Treatment, Year) %>% 
  filter(Seeded == "Summer" & Treatment != "NON") %>% 
  summarise(max_cover = sum(max_cover))
## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.
SR_summer_zero <- SR_summer_zero  %>% 
  filter(max_cover == 0) 


SR_summer_zero <- SR_summer_zero [,c(1,2,3)]
SR_summer_zero$Summer_SR <- rep(0, nrow(SR_summer_zero))


SR_summer_all <- full_join(SR_summer,SR_summer_zero)
## Joining with `by = join_by(Block, Treatment, Year, Summer_SR)`
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Summer_SR ~ Treatment + Year + (1 | Block)
##    Data: SR_summer_all
## 
## REML criterion at convergence: 133.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3183 -0.5433 -0.2825  0.5178  2.1606 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 0.06321  0.2514  
##  Residual             0.92763  0.9631  
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    2.2917     0.3274 34.2203   7.000 4.31e-08 ***
## TreatmentLL   -1.9167     0.3932 38.0000  -4.875 1.96e-05 ***
## TreatmentNAT  -1.9167     0.3932 38.0000  -4.875 1.96e-05 ***
## TreatmentSIM  -0.1667     0.3932 38.0000  -0.424    0.674    
## Year2023       0.2500     0.2780 38.0000   0.899    0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.601                     
## TreatmntNAT -0.601  0.500              
## TreatmntSIM -0.601  0.500  0.500       
## Year2023    -0.425  0.000  0.000  0.000
##            R2m       R2c
## [1,] 0.4697125 0.5035431
Anova(SR.summer.mod.normal, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Summer_SR
##               Chisq Df Pr(>Chisq)    
## (Intercept) 49.0067  1  2.551e-12 ***
## Treatment   43.6596  3  1.783e-09 ***
## Year         0.8085  1     0.3686    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SR.summer.mod.normal.est <- lmer(Summer_SR~Treatment+(1|Block), data =
SR_summer_all)


est.lmer.SR.summer.treatment <- emmeans::emmeans(SR.summer.mod.normal.est, ~ Treatment, type = "response")

pairs(est.lmer.SR.summer.treatment , adjust = "tukey")
##  contrast  estimate    SE df t.ratio p.value
##  LE - LL      1.917 0.392 39   4.887  0.0001
##  LE - NAT     1.917 0.392 39   4.887  0.0001
##  LE - SIM     0.167 0.392 39   0.425  0.9739
##  LL - NAT     0.000 0.392 39   0.000  1.0000
##  LL - SIM    -1.750 0.392 39  -4.462  0.0004
##  NAT - SIM   -1.750 0.392 39  -4.462  0.0004
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates
set.seed(7)

SR_summer_all_mean_sd <- SR_summer_all %>% 
  group_by(Treatment, Year) %>% 
  summarize(mean_SR_summer = mean(Summer_SR), sd_SR_summer = sd(Summer_SR))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
SR_summer_all_mean_sd$Treatment <- factor(SR_summer_all_mean_sd$Treatment, levels=c("NAT", "LE", "LL", "SIM"))

SR_summer_all$Treatment <- factor(SR_summer_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))

sr.summer.plot_treatment <-   ggplot()+
  geom_jitter(data = SR_summer_all, 
              aes(x = Year, y = Summer_SR, fill = Treatment, 
                  shape = Treatment, group = Treatment), 
              alpha = 1, size = 3.5, 
              position = position_jitterdodge(dodge.width = 0.85))+
  geom_errorbar(data = SR_summer_all_mean_sd,
                aes(x = Year, y = mean_SR_summer,  group = Treatment, 
                    ymin = mean_SR_summer - sd_SR_summer, 
                    ymax = mean_SR_summer + sd_SR_summer), 
                position = position_dodge(width = 0.85), width = .5,
                size = 0.75)+
  geom_point(data = SR_summer_all_mean_sd, size = 4.5, 
             aes(x = Year, y = mean_SR_summer, 
                 shape = Treatment,  group = Treatment), 
             position = position_dodge2(width = 0.85),
             fill = "black")+
  labs(y = "Sown summer richness", x = "" )+
  theme_classic() +
  theme(text=element_text(size=18), 
        legend.key.size=unit(1, "cm"),
        legend.position="none",  
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 16))+
  scale_fill_manual(name = "Treatment",
                    breaks = c("LE", "LL", "NAT", "SIM"),  
                    labels = c("Summer first", "Fall first", "Natural", "Simultaneous"), 
                    values = c("#009E73", "#56B4E9",  "#F0E442", "#0072B2"))+
  # scale_x_discrete(breaks = c("LE", "LL", "NAT", "SIM"),
  #                  labels = c("Summer first", "Fall first", "Natural", "Simultaneous"))+
  scale_shape_manual(name = "Treatment", 
                     breaks = c("LE", "LL", "NAT", "SIM"), 
                     labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
                     values=c(21, 22,23, 25))+
  scale_y_continuous(breaks = seq(0, 20, by = 5),
                     limits = c(-0.75, 20))




sr.summer.plot_treatment <- sr.summer.plot_treatment +
  annotate("text", x = 1.2, y = 20, label = "Treatment: p < 0.001***",
           fontface = 2, size = 4)

Fall-Dispersing

### Filter out everything that isn't a seeded fall-dispersing species

SR_fall <- Sown_species_ID  %>% 
  filter(Seeded == "Fall" & max_cover > 0 & Treatment != "NON") %>% 
  summarize(Fall_SR = n_distinct(SPP6))
## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.
SR_fall_zero <- Sown_species_ID %>% 
  group_by(Block, Treatment, Year) %>% 
  filter(Seeded == "Fall"  & Treatment != "NON") %>% 
  summarise(max_cover = sum(max_cover))
## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.
SR_fall_zero <- SR_fall_zero  %>% 
  filter(max_cover == 0) 


SR_fall_zero <- SR_fall_zero[,c(1,2,3)]
SR_fall_zero$Fall_SR <- rep(0, nrow(SR_fall_zero))


SR_fall_all <- full_join(SR_fall,SR_fall_zero)
## Joining with `by = join_by(Block, Treatment, Year, Fall_SR)`
## 
## Call:
## lm(formula = Fall_SR ~ Treatment + Year, data = SR_fall_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9583 -0.7083  0.0417  0.7292  3.3750 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.7083     0.4977   5.442 2.36e-06 ***
## TreatmentLL    3.8333     0.6295   6.089 2.71e-07 ***
## TreatmentNAT   0.5000     0.6295   0.794    0.431    
## TreatmentSIM   5.1667     0.6295   8.207 2.42e-10 ***
## Year2023      -0.2500     0.4452  -0.562    0.577    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.542 on 43 degrees of freedom
## Multiple R-squared:  0.6919, Adjusted R-squared:  0.6633 
## F-statistic: 24.15 on 4 and 43 DF,  p-value: 1.608e-10
Anova(SR.fall.mod.normal , type = 3)
## Anova Table (Type III tests)
## 
## Response: Fall_SR
##              Sum Sq Df F value    Pr(>F)    
## (Intercept)  70.417  1 29.6129 2.356e-06 ***
## Treatment   228.917  3 32.0894 4.753e-11 ***
## Year          0.750  1  0.3154    0.5773    
## Residuals   102.250 43                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SR.fall.mod.normal.est <- lm(Fall_SR~Treatment, data = SR_fall_all)


est.lmer.SR.fall.treatment <- emmeans::emmeans(SR.fall.mod.normal.est, ~ Treatment, type = "response")

pairs(est.lmer.SR.fall.treatment , adjust = "tukey")
##  contrast  estimate    SE df t.ratio p.value
##  LE - LL      -3.83 0.625 44  -6.137  <.0001
##  LE - NAT     -0.50 0.625 44  -0.800  0.8538
##  LE - SIM     -5.17 0.625 44  -8.272  <.0001
##  LL - NAT      3.33 0.625 44   5.337  <.0001
##  LL - SIM     -1.33 0.625 44  -2.135  0.1582
##  NAT - SIM    -4.67 0.625 44  -7.471  <.0001
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
set.seed(5)

SR_fall_all_mean_sd <- SR_fall_all %>% 
  group_by(Treatment, Year) %>% 
  summarize(mean_SR_fall = mean(Fall_SR), sd_SR_fall = sd(Fall_SR))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
SR_fall_all_mean_sd $Treatment <- factor(SR_fall_all_mean_sd $Treatment, levels=c("NAT", "LE", "LL", "SIM"))

SR_fall_all$Treatment <- factor(SR_fall_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))

sr.fall.plot_treatment <-   ggplot()+
    geom_jitter(data = SR_fall_all,
              aes(x = Year, y = Fall_SR, fill = Treatment, shape = Treatment, group = Treatment), 
              alpha = 0.8, size = 3.5,
              position = position_jitterdodge(dodge.width = 0.85))+
    geom_errorbar(data = SR_fall_all_mean_sd, 
                  aes(x = Year, y = mean_SR_fall, ymin = mean_SR_fall - sd_SR_fall, 
                      ymax = mean_SR_fall + sd_SR_fall, group = Treatment),
                  position = position_dodge(width = 0.85), width = .5, size = 0.75) +
    geom_point(data = SR_fall_all_mean_sd, size = 4.5, 
               aes(x = Year, y = mean_SR_fall, shape = Treatment, group = Treatment), 
               fill = "black", position = position_dodge(width = 0.85))+
    labs(y = "Sown fall richness", x = "" ) +
    theme_classic() +
    theme(text=element_text(size=18), 
          legend.key.size=unit(1, "cm"),
          legend.position="none",  
          axis.text.x = element_text(size = 16), 
          axis.text.y = element_text(size = 16))+
    scale_fill_manual(name = "Treatment", 
                      breaks = c("LE", "LL", "NAT", "SIM"), 
                      labels = c("Summer first", "Fall first", "Natural", "Simultaneous"), 
                      values = c("#009E73", "#56B4E9",  "#F0E442", "#0072B2"))+
    scale_shape_manual(name = "Treatment", 
                       breaks = c("LE", "LL", "NAT", "SIM"),
                       labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
                       values=c(21, 22,23, 25))+
  scale_y_continuous(breaks = seq(0, 20, by = 5),
                     limits = c(-0.75, 20))



sr.fall.plot_treatment  <- sr.fall.plot_treatment +
  annotate("text", x = 1.2, y = 20, label = "Treatment: p < 0.001***",
           fontface = 2, size = 4)

Unsown

### Filter out everything that isn't a seeded fall-dispersing species

SR_unsown_all <- Sown_species_ID  %>% 
  filter(Seeded == "Unseeded" & max_cover > 0 & Treatment != "NON") %>% 
  summarize(unsown_SR = n_distinct(SPP6))
## `summarise()` has grouped output by 'Block', 'Treatment'. You can override
## using the `.groups` argument.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: unsown_SR ~ Treatment + Year + (1 | Block)
##    Data: SR_unsown_all
## 
## REML criterion at convergence: 212.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.11194 -0.60308  0.06422  0.64409  1.79134 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 0.5141   0.717   
##  Residual             5.8202   2.413   
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   15.0833     0.8318 32.4490  18.132  < 2e-16 ***
## TreatmentLL   -1.5833     0.9849 38.0000  -1.608  0.11620    
## TreatmentNAT  -0.3333     0.9849 38.0000  -0.338  0.73689    
## TreatmentSIM  -0.7500     0.9849 38.0000  -0.761  0.45106    
## Year2023      -2.0000     0.6964 38.0000  -2.872  0.00664 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.592                     
## TreatmntNAT -0.592  0.500              
## TreatmntSIM -0.592  0.500  0.500       
## Year2023    -0.419  0.000  0.000  0.000
Anova(SR.unsown.mod.normal , type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: unsown_SR
##                Chisq Df Pr(>Chisq)    
## (Intercept) 328.7867  1  < 2.2e-16 ***
## Treatment     2.8922  3   0.408540    
## Year          8.2472  1   0.004082 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
set.seed(5)

SR_unsown_all_mean_sd <- SR_unsown_all %>% 
  group_by(Treatment, Year) %>% 
  summarize(mean_SR_unsown = mean(unsown_SR), sd_SR_unsown = sd(unsown_SR))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
SR_unsown_all_mean_sd $Treatment <- factor(SR_unsown_all_mean_sd $Treatment, levels=c("NAT", "LE", "LL", "SIM"))

SR_unsown_all$Treatment <- factor(SR_unsown_all$Treatment, levels=c("NAT", "LE", "LL", "SIM"))

sr.unsown.plot_treatment <-   ggplot()+
    geom_jitter(data = SR_unsown_all,
              aes(x = Year, y = unsown_SR, fill = Treatment, shape = Treatment, group = Treatment), 
              alpha = 0.8, size = 3.5,
              position = position_jitterdodge(dodge.width = 0.85))+
    geom_errorbar(data = SR_unsown_all_mean_sd, 
                  aes(x = Year, y = mean_SR_unsown, ymin = mean_SR_unsown - sd_SR_unsown, 
                      ymax = mean_SR_unsown + sd_SR_unsown, group = Treatment),
                  position = position_dodge(width = 0.85), width = .5, size = 0.75) +
    geom_point(data = SR_unsown_all_mean_sd, size = 4.5, 
               aes(x = Year, y = mean_SR_unsown, shape = Treatment, group = Treatment), 
               fill = "black", position = position_dodge(width = 0.85))+
    labs(y = "Unsown richness", x = "" ) +
    theme_classic() +
    theme(text=element_text(size=18), 
          legend.key.size=unit(1, "cm"),
          legend.position="none",  
          axis.text.x = element_text(size = 16), 
          axis.text.y = element_text(size = 16))+
    scale_fill_manual(name = "Treatment", 
                      breaks = c("LE", "LL", "NAT", "SIM"), 
                      labels = c("Summer first", "Fall first", "Natural", "Simultaneous"), 
                      values = c("#009E73", "#56B4E9",  "#F0E442", "#0072B2"))+
    scale_shape_manual(name = "Treatment", 
                       breaks = c("LE", "LL", "NAT", "SIM"),
                       labels = c("Summer first", "Fall first", "Natural", "Simultaneous"),
                       values=c(21, 22,23, 25))+
  scale_y_continuous(breaks = seq(0, 20, by = 5),
                     limits = c(-0.75, 20))


sr.unsown.plot_treatment

Shannon Diversity Index

SDI - Model assumptions

SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.

SDI - Model fitting

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
##    Data: Summary_table
## 
## REML criterion at convergence: 7.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.42601 -0.49701  0.08675  0.60144  1.94713 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block    (Intercept) 0.004618 0.06796 
##  Residual             0.049227 0.22187 
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.026960   0.076796 31.983342  26.394  < 2e-16 ***
## TreatmentLE   0.122187   0.090579 38.000000   1.349  0.18533    
## TreatmentLL   0.173774   0.090579 38.000000   1.918  0.06259 .  
## TreatmentSIM  0.251361   0.090579 38.000000   2.775  0.00851 ** 
## Year2023      0.008802   0.064049 38.000000   0.137  0.89142    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLE TrtmLL TrtSIM
## TreatmentLE -0.590                     
## TreatmentLL -0.590  0.500              
## TreatmntSIM -0.590  0.500  0.500       
## Year2023    -0.417  0.000  0.000  0.000
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Shannon
##                Chisq Df Pr(>Chisq)    
## (Intercept) 696.6555  1    < 2e-16 ***
## Treatment     8.1465  3    0.04308 *  
## Year          0.0189  1    0.89070    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##            R2m       R2c
## [1,] 0.1370613 0.2110774
##  contrast  estimate     SE df t.ratio p.value
##  NAT - LE   -0.1222 0.0906 38  -1.349  0.5384
##  NAT - LL   -0.1738 0.0906 38  -1.918  0.2375
##  NAT - SIM  -0.2514 0.0906 38  -2.775  0.0406
##  LE - LL    -0.0516 0.0906 38  -0.570  0.9405
##  LE - SIM   -0.1292 0.0906 38  -1.426  0.4913
##  LL - SIM   -0.0776 0.0906 38  -0.857  0.8269
## 
## Results are averaged over the levels of: Year 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

SDI - Figures

Plot Panel

Question 2: Does timing of species arrival influence seeded cover?

note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)

Seeded2 <- Seeded %>% filter(Year == "2023") %>% filter(max_cover > 0) 
unique(Seeded2$SPP6)
##  [1] "ACHMIL" "BIDARI" "CHAFAS" "CORLAN" "CORTRI" "ERYYUC" "HELMOL" "KOEMAC"
##  [9] "LESCAP" "LIAPYC" "MONFIS" "PENDIG" "PYCTEN" "RATPIN" "RUDHIR" "SORNUT"
## [17] "SPHOBT" "SPOHET" "VIOPED"

Total Seeded Cover

TSC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year, 
##     family = binomial, data = Total_Cover_Tot)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -2.11998    0.08881 -23.870  < 2e-16 ***
## TreatmentLL            0.94324    0.10412   9.059  < 2e-16 ***
## TreatmentNAT          -0.38112    0.13799  -2.762  0.00575 ** 
## TreatmentSIM           1.06270    0.10397  10.221  < 2e-16 ***
## Year2023               0.62783    0.11282   5.565 2.63e-08 ***
## TreatmentLL:Year2023  -0.41302    0.13685  -3.018  0.00254 ** 
## TreatmentNAT:Year2023 -0.36533    0.18174  -2.010  0.04441 *  
## TreatmentSIM:Year2023 -0.37777    0.13662  -2.765  0.00569 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 837.79  on 47  degrees of freedom
## Residual deviance: 293.23  on 40  degrees of freedom
## AIC: 557.58
## 
## Number of Fisher Scoring iterations: 5
## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                LR Chisq Df Pr(>Chisq)    
## Treatment       272.445  3  < 2.2e-16 ***
## Year             32.035  1  1.514e-08 ***
## Treatment:Year   10.305  3    0.01614 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast                    odds.ratio     SE  df null z.ratio p.value
##  LE Year2022 / LL Year2022        0.389 0.0405 Inf    1  -9.059  <.0001
##  LE Year2022 / NAT Year2022       1.464 0.2020 Inf    1   2.762  0.1048
##  LE Year2022 / SIM Year2022       0.346 0.0359 Inf    1 -10.221  <.0001
##  LE Year2022 / LE Year2023        0.534 0.0602 Inf    1  -5.565  <.0001
##  LE Year2022 / LL Year2023        0.314 0.0328 Inf    1 -11.076  <.0001
##  LE Year2022 / NAT Year2023       1.126 0.1469 Inf    1   0.909  0.9853
##  LE Year2022 / SIM Year2023       0.269 0.0281 Inf    1 -12.574  <.0001
##  LL Year2022 / NAT Year2022       3.760 0.4466 Inf    1  11.150  <.0001
##  LL Year2022 / SIM Year2022       0.887 0.0680 Inf    1  -1.558  0.7751
##  LL Year2022 / LE Year2023        1.371 0.1210 Inf    1   3.572  0.0085
##  LL Year2022 / LL Year2023        0.807 0.0625 Inf    1  -2.774  0.1016
##  LL Year2022 / NAT Year2023       2.892 0.3181 Inf    1   9.654  <.0001
##  LL Year2022 / SIM Year2023       0.691 0.0534 Inf    1  -4.784  <.0001
##  NAT Year2022 / SIM Year2022      0.236 0.0280 Inf    1 -12.169  <.0001
##  NAT Year2022 / LE Year2023       0.365 0.0461 Inf    1  -7.977  <.0001
##  NAT Year2022 / LL Year2023       0.215 0.0256 Inf    1 -12.917  <.0001
##  NAT Year2022 / NAT Year2023      0.769 0.1096 Inf    1  -1.842  0.5909
##  NAT Year2022 / SIM Year2023      0.184 0.0219 Inf    1 -14.231  <.0001
##  SIM Year2022 / LE Year2023       1.545 0.1361 Inf    1   4.935  <.0001
##  SIM Year2022 / LL Year2023       0.909 0.0702 Inf    1  -1.234  0.9218
##  SIM Year2022 / NAT Year2023      3.259 0.3580 Inf    1  10.753  <.0001
##  SIM Year2022 / SIM Year2023      0.779 0.0600 Inf    1  -3.246  0.0258
##  LE Year2023 / LL Year2023        0.588 0.0523 Inf    1  -5.971  <.0001
##  LE Year2023 / NAT Year2023       2.109 0.2495 Inf    1   6.312  <.0001
##  LE Year2023 / SIM Year2023       0.504 0.0447 Inf    1  -7.729  <.0001
##  LL Year2023 / NAT Year2023       3.585 0.3958 Inf    1  11.563  <.0001
##  LL Year2023 / SIM Year2023       0.857 0.0667 Inf    1  -1.988  0.4901
##  NAT Year2023 / SIM Year2023      0.239 0.0264 Inf    1 -12.981  <.0001
## 
## P value adjustment: tukey method for comparing a family of 8 estimates 
## Tests are performed on the log odds ratio scale
##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL        0.479 0.0328 Inf    1 -10.767  <.0001
##  LE / NAT       1.757 0.1597 Inf    1   6.204  <.0001
##  LE / SIM       0.417 0.0285 Inf    1 -12.792  <.0001
##  LL / NAT       3.671 0.2977 Inf    1  16.039  <.0001
##  LL / SIM       0.872 0.0476 Inf    1  -2.510  0.0584
##  NAT / SIM      0.237 0.0192 Inf    1 -17.751  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale
##  contrast            odds.ratio     SE  df null z.ratio p.value
##  Year2022 / Year2023      0.713 0.0378 Inf    1  -6.391  <.0001
## 
## Results are averaged over the levels of: Treatment 
## Tests are performed on the log odds ratio scale

TSC - Figures

Seeded cover of summer dispersing species

ESC - Model fitting

## 
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year, 
##     family = binomial, data = Total_Cover_Early)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.9526     0.1056 -27.955  < 2e-16 ***
## TreatmentLL   -3.3561     0.3411  -9.838  < 2e-16 ***
## TreatmentNAT  -2.0984     0.1995 -10.519  < 2e-16 ***
## TreatmentSIM  -0.3623     0.1070  -3.387 0.000707 ***
## Year2023       1.0821     0.1106   9.780  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 734.86  on 47  degrees of freedom
## Residual deviance: 247.79  on 43  degrees of freedom
## AIC: 381.05
## 
## Number of Fisher Scoring iterations: 6
## Analysis of Deviance Table (Type III tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##           LR Chisq Df Pr(>Chisq)    
## Treatment   373.50  3  < 2.2e-16 ***
## Year        106.69  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL      28.6759 9.7820 Inf    1   9.838  <.0001
##  LE / NAT      8.1535 1.6265 Inf    1  10.519  <.0001
##  LE / SIM      1.4366 0.1537 Inf    1   3.387  0.0039
##  LL / NAT      0.2843 0.1089 Inf    1  -3.284  0.0056
##  LL / SIM      0.0501 0.0172 Inf    1  -8.706  <.0001
##  NAT / SIM     0.1762 0.0360 Inf    1  -8.502  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale
##  contrast            odds.ratio     SE  df null z.ratio p.value
##  Year2022 / Year2023      0.339 0.0375 Inf    1  -9.780  <.0001
## 
## Results are averaged over the levels of: Treatment 
## Tests are performed on the log odds ratio scale

ESC - Figures

Fall dispersing species only

LSC - Model fitting

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
##    Data: Total_Cover_Late
## 
##      AIC      BIC   logLik deviance df.resid 
##    605.3    616.5   -296.6    593.3       42 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3147 -1.7714 -0.6057  1.2574  6.5425 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block  (Intercept) 0.008129 0.09016 
## Number of obs: 48, groups:  Block, 6
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.74863    0.09364 -29.352  < 2e-16 ***
## TreatmentLL   1.59358    0.09078  17.554  < 2e-16 ***
## TreatmentNAT  0.16919    0.11161   1.516  0.12953    
## TreatmentSIM  1.55200    0.09191  16.885  < 2e-16 ***
## Year2023      0.13908    0.05056   2.751  0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796                     
## TreatmntNAT -0.643  0.664              
## TreatmntSIM -0.789  0.807  0.656       
## Year2023    -0.279  0.018 -0.002  0.024
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: cbind(seeded_cover, tot_cover)
##                Chisq Df Pr(>Chisq)    
## (Intercept) 861.5598  1    < 2e-16 ***
## Treatment   566.5093  3    < 2e-16 ***
## Year          7.5683  1    0.00594 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  contrast  odds.ratio     SE  df null z.ratio p.value
##  LE / LL        0.203 0.0184 Inf    1 -17.554  <.0001
##  LE / NAT       0.844 0.0942 Inf    1  -1.516  0.4278
##  LE / SIM       0.212 0.0195 Inf    1 -16.885  <.0001
##  LL / NAT       4.155 0.3536 Inf    1  16.738  <.0001
##  LL / SIM       1.042 0.0592 Inf    1   0.732  0.8841
##  NAT / SIM      0.251 0.0216 Inf    1 -16.031  <.0001
## 
## Results are averaged over the levels of: Year 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale
##  contrast            odds.ratio    SE  df null z.ratio p.value
##  Year2022 / Year2023       0.87 0.044 Inf    1  -2.751  0.0059
## 
## Results are averaged over the levels of: Treatment 
## Tests are performed on the log odds ratio scale

LSC - Figures

Plot Panel

Species specific differences

## 
##  One Sample t-test
## 
## data:  test2$difference
## t = 5.4921, df = 89, p-value = 3.722e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3.935637 8.397697
## sample estimates:
## mean of x 
##  6.166667

Question 3: Does timing of species arrival influence species composition? - Bray-Curtis

Year: 2022

## Call: capscale(formula = inner.hel.2022 ~ Treatment, data =
## env.data.2022, distance = "bray", add = TRUE)
## 
##               Inertia Proportion Rank
## Total          3.7536     1.0000     
## Constrained    0.8657     0.2306    3
## Unconstrained  2.8879     0.7694   20
## Inertia is Lingoes adjusted squared Bray distance 
## Species scores projected from 'inner.hel.2022' 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3 
## 0.5388 0.2425 0.0843 
## 
## Eigenvalues for unconstrained axes:
##   MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8 
## 0.5403 0.4052 0.2960 0.2529 0.1945 0.1649 0.1568 0.1338 
## (Showing 8 of 20 unconstrained eigenvalues)
## 
## Constant added to distances: 0.03513658
## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = inner.hel.2022 ~ Treatment, data = env.data.2022, distance = "bray", add = TRUE)
##           Df SumOfSqs      F Pr(>F)    
## Treatment  3  0.86567 1.9984  0.001 ***
## Residual  20  2.88789                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Multiple comparisons for capscale for all contrasts of Treatment 
## 
## Model:  BiodiversityR::multiconstrained(method = "capscale", formula = inner.hel.2022 ~ Treatment, data = env.data.2022, distance = "bray", add = TRUE, by = "term") 
## 
##             Df SumOfSqs      F Pr(>F)   
## LE vs. LL    1  0.32432 3.2161  0.002 **
## LE vs. NAT   1  0.08466 0.6738  0.812   
## LE vs. SIM   1  0.31377 2.6708  0.007 **
## LL vs. NAT   1  0.25386 2.4115  0.006 **
## LL vs. SIM   1  0.17717 1.9078  0.030 * 
## NAT vs. SIM  1  0.37810 3.0297  0.006 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   UniqueBlock Block Treatment Year Bare_cover  Above Ground Middle       axis1
## 1        LE_1     1        LE <NA>         NA 1575.5     67  239.5 -0.07771696
## 2        LL_1     1        LL 2022         14 1565.0     31  215.0  0.48336017
## 3       NAT_1     1       NAT 2022          4 1535.0     32  128.5 -0.43746584
## 4       SIM_1     1       SIM 2022          2 1580.5     22  626.0  0.58299824
## 5        LE_2     2        LE 2022         10 1611.0    163  422.0 -0.61581200
## 6        LL_2     2        LL 2022         11 1601.5     71  336.0  0.27390902
##         axis2 labels
## 1  0.38755734   LE_1
## 2 -0.95497148   LL_1
## 3 -0.03453055  NAT_1
## 4  1.97668133  SIM_1
## 5  1.09725065   LE_2
## 6 -0.60130836   LL_2
##                   axis1         axis2     labels
## ACAVIR     -0.013759760  0.0105436448     ACAVIR
## ACHMIL     -0.039905108  0.3151005315     ACHMIL
## AGRHYE      0.008944552  0.0420245595     AGRHYE
## AMBART     -0.165234521  0.0893240283     AMBART
## AMBTRI     -0.395292407 -0.4196119068     AMBTRI
## BARVUL     -0.022532498 -0.0940021827     BARVUL
## BIDARI      0.616701687 -0.0381609416     BIDARI
## BROJAP      0.113385176  0.0873793317     BROJAP
## CARSPP     -0.133319955 -0.0709646462     CARSPP
## CERSPP      0.013460740  0.0070992909     CERSPP
## CHAFAS      0.221973715 -0.0615936378     CHAFAS
## CONCAN      0.026314578  0.0186473255     CONCAN
## CORLAN      0.047431342  0.1387802756     CORLAN
## CORPAL     -0.043650430  0.0230115005     CORPAL
## CORTRI      0.259290736  0.0956570002     CORTRI
## DESSPP     -0.002361261 -0.0276468238     DESSPP
## DIGISC     -0.005805795  0.0523668346     DIGISC
## DIGSAN      0.004596870  0.0292496220     DIGSAN
## ECHCRU     -0.007543745 -0.0139571869     ECHCRU
## ERASPE      0.017765474  0.0161270608     ERASPE
## ERIANN     -0.020170112 -0.0512415405     ERIANN
## ERYYUC      0.034816004 -0.0109188717     ERYYUC
## FESARU      0.001119958 -0.0145731064     FESARU
## GALPED     -0.008890104  0.0004354018     GALPED
## GERCAR     -0.029033113 -0.0210649548     GERCAR
## GEUCAN     -0.018517585 -0.0096680542     GEUCAN
## HELAUT     -0.105661417 -0.0149522445     HELAUT
## HELMOL      0.129128316 -0.0683750610     HELMOL
## HORPUS      0.035115845 -0.0733476653     HORPUS
## IPOLAC     -0.006344513  0.0278887241     IPOLAC
## IVAANN      0.262803154 -0.0578486576     IVAANN
## JUNINT      0.011148351 -0.0176837652     JUNINT
## LACSER      0.017301947  0.0157062823     LACSER
## LESCAP      0.116802974 -0.0116205743     LESCAP
## LIAPYC      0.017301947  0.0157062823     LIAPYC
## LINSUL     -0.013759760  0.0105436448     LINSUL
## LONMAC      0.004421824  0.0270484397     LONMAC
## LOTCOR     -0.082750521  0.3856035839     LOTCOR
## MELSPP     -0.043445692 -0.0226830500     MELSPP
## MONFIS     -0.009543532  0.0202951310     MONFIS
## MYOVER     -0.039536260  0.0611374149     MYOVER
## OXADIL     -0.032682355  0.0243796838     OXADIL
## PANCAP      0.003201855  0.0279832615     PANCAP
## PANDIC     -0.003845438 -0.0200239744     PANDIC
## PENDIG     -0.082854822  0.0077086447     PENDIG
## PERHYD/PUN -0.015929940  0.0122065818 PERHYD/PUN
## PHAARU     -0.106205632 -0.0554500935     PHAARU
## PLALAN      0.037094780 -0.0381481848     PLALAN
## PLAMAJ     -0.037542985  0.0089781276     PLAMAJ
## PLAVIR     -0.036194605 -0.0025266870     PLAVIR
## POAPRA      0.046096825  0.0669222606     POAPRA
## POTNOR      0.034618018 -0.0117608580     POTNOR
## RANABO      0.063310648  0.0623983669     RANABO
## RATPIN      0.295359023 -0.1903300280     RATPIN
## RUDHIR      0.128747186 -0.1095167813     RUDHIR
## RUMCRI     -0.128181852  0.0250011764     RUMCRI
## SETFAB     -0.087935392 -0.0525684638     SETFAB
## SETPUM     -0.027218251  0.0588510004     SETPUM
## SIDSPI      0.019576599  0.0177711555     SIDSPI
## SOLCAN     -0.015360372 -0.0080196692     SOLCAN
## SOLCAR      0.013348165 -0.0229292736     SOLCAR
## SORNUT      0.423351229 -0.0195015817     SORNUT
## SPHOBT      0.031491428  0.0285871434     SPHOBT
## SPOHET      0.025712643  0.0233413048     SPOHET
## TAROFF      0.005036684  0.0524681788     TAROFF
## TRIPRA      0.040655258  0.0369058434     TRIPRA
## TRIREP     -0.379489398 -0.0327028056     TRIREP
## VERARV     -0.035216468  0.0164490440     VERARV
## VERURT      0.019309511 -0.0306291798     VERURT
## VIOPED     -0.017935072  0.0434831986     VIOPED
##   axis     ggplot        label
## 1    1 xlab.label CAP1 (14.4%)
## 2    2 ylab.label  CAP2 (6.5%)
##                     r     p        axis1         axis2     labels
## ACAVIR     0.01156079 0.916 -0.013759760  0.0105436448     ACAVIR
## ACHMIL     0.26416995 0.036 -0.039905108  0.3151005315     ACHMIL
## AGRHYE     0.04221726 0.640  0.008944552  0.0420245595     AGRHYE
## AMBART     0.05133255 0.574 -0.165234521  0.0893240283     AMBART
## AMBTRI     0.62167784 0.001 -0.395292407 -0.4196119068     AMBTRI
## BARVUL     0.24118220 0.062 -0.022532498 -0.0940021827     BARVUL
## BIDARI     0.80347542 0.001  0.616701687 -0.0381609416     BIDARI
## BROJAP     0.11728601 0.306  0.113385176  0.0873793317     BROJAP
## CARSPP     0.08699372 0.381 -0.133319955 -0.0709646462     CARSPP
## CERSPP     0.01879381 0.846  0.013460740  0.0070992909     CERSPP
## CHAFAS     0.32314850 0.016  0.221973715 -0.0615936378     CHAFAS
## CONCAN     0.05347984 0.563  0.026314578  0.0186473255     CONCAN
## CORLAN     0.23679300 0.068  0.047431342  0.1387802756     CORLAN
## CORPAL     0.22152803 0.072 -0.043650430  0.0230115005     CORPAL
## CORTRI     0.57092202 0.002  0.259290736  0.0956570002     CORTRI
## DESSPP     0.08205552 0.426 -0.002361261 -0.0276468238     DESSPP
## DIGISC     0.06806011 0.495 -0.005805795  0.0523668346     DIGISC
## DIGSAN     0.31708125 0.016  0.004596870  0.0292496220     DIGSAN
## ECHCRU     0.01043125 0.912 -0.007543745 -0.0139571869     ECHCRU
## ERASPE     0.16914332 0.090  0.017765474  0.0161270608     ERASPE
## ERIANN     0.10723478 0.337 -0.020170112 -0.0512415405     ERIANN
## ERYYUC     0.17864825 0.129  0.034816004 -0.0109188717     ERYYUC
## FESARU     0.01976619 0.825  0.001119958 -0.0145731064     FESARU
## GALPED     0.07257478 0.453 -0.008890104  0.0004354018     GALPED
## GERCAR     0.03530864 0.695 -0.029033113 -0.0210649548     GERCAR
## GEUCAN     0.07878386 0.578 -0.018517585 -0.0096680542     GEUCAN
## HELAUT     0.02234030 0.824 -0.105661417 -0.0149522445     HELAUT
## HELMOL     0.32842446 0.010  0.129128316 -0.0683750610     HELMOL
## HORPUS     0.09919497 0.349  0.035115845 -0.0733476653     HORPUS
## IPOLAC     0.09105548 0.401 -0.006344513  0.0278887241     IPOLAC
## IVAANN     0.44187949 0.004  0.262803154 -0.0578486576     IVAANN
## JUNINT     0.08948316 0.472  0.011148351 -0.0176837652     JUNINT
## LACSER     0.01396845 0.859  0.017301947  0.0157062823     LACSER
## LESCAP     0.24964718 0.045  0.116802974 -0.0116205743     LESCAP
## LIAPYC     0.01396845 0.859  0.017301947  0.0157062823     LIAPYC
## LINSUL     0.01156079 0.916 -0.013759760  0.0105436448     LINSUL
## LONMAC     0.09228958 0.362  0.004421824  0.0270484397     LONMAC
## LOTCOR     0.60180742 0.001 -0.082750521  0.3856035839     LOTCOR
## MELSPP     0.12993356 0.128 -0.043445692 -0.0226830500     MELSPP
## MONFIS     0.09716618 0.320 -0.009543532  0.0202951310     MONFIS
## MYOVER     0.10341326 0.323 -0.039536260  0.0611374149     MYOVER
## OXADIL     0.04738271 0.598 -0.032682355  0.0243796838     OXADIL
## PANCAP     0.13538690 0.231  0.003201855  0.0279832615     PANCAP
## PANDIC     0.04478159 0.623 -0.003845438 -0.0200239744     PANDIC
## PENDIG     0.09639392 0.349 -0.082854822  0.0077086447     PENDIG
## PERHYD/PUN 0.05224616 0.742 -0.015929940  0.0122065818 PERHYD/PUN
## PHAARU     0.01682139 0.855 -0.106205632 -0.0554500935     PHAARU
## PLALAN     0.07651040 0.445  0.037094780 -0.0381481848     PLALAN
## PLAMAJ     0.13775527 0.204 -0.037542985  0.0089781276     PLAMAJ
## PLAVIR     0.03838884 0.702 -0.036194605 -0.0025266870     PLAVIR
## POAPRA     0.36792330 0.006  0.046096825  0.0669222606     POAPRA
## POTNOR     0.02307529 0.791  0.034618018 -0.0117608580     POTNOR
## RANABO     0.36519063 0.007  0.063310648  0.0623983669     RANABO
## RATPIN     0.75052506 0.001  0.295359023 -0.1903300280     RATPIN
## RUDHIR     0.22610230 0.069  0.128747186 -0.1095167813     RUDHIR
## RUMCRI     0.14928772 0.185 -0.128181852  0.0250011764     RUMCRI
## SETFAB     0.06697542 0.477 -0.087935392 -0.0525684638     SETFAB
## SETPUM     0.08685742 0.384 -0.027218251  0.0588510004     SETPUM
## SIDSPI     0.30323356 0.048  0.019576599  0.0177711555     SIDSPI
## SOLCAN     0.12993356 0.128 -0.015360372 -0.0080196692     SOLCAN
## SOLCAR     0.02744986 0.774  0.013348165 -0.0229292736     SOLCAR
## SORNUT     0.70787270 0.001  0.423351229 -0.0195015817     SORNUT
## SPHOBT     0.11670926 0.333  0.031491428  0.0285871434     SPHOBT
## SPOHET     0.11670926 0.333  0.025712643  0.0233413048     SPOHET
## TAROFF     0.08199427 0.406  0.005036684  0.0524681788     TAROFF
## TRIPRA     0.11670926 0.333  0.040655258  0.0369058434     TRIPRA
## TRIREP     0.45984733 0.002 -0.379489398 -0.0327028056     TRIREP
## VERARV     0.03991363 0.650 -0.035216468  0.0164490440     VERARV
## VERURT     0.08948316 0.472  0.019309511 -0.0306291798     VERURT
## VIOPED     0.01282063 0.878 -0.017935072  0.0434831986     VIOPED
##                r     p       axis1       axis2 labels
## AMBTRI 0.6216778 0.001 -0.39529241 -0.41961191 AMBTRI
## BIDARI 0.8034754 0.001  0.61670169 -0.03816094 BIDARI
## CORTRI 0.5709220 0.002  0.25929074  0.09565700 CORTRI
## IVAANN 0.4418795 0.004  0.26280315 -0.05784866 IVAANN
## LOTCOR 0.6018074 0.001 -0.08275052  0.38560358 LOTCOR
## POAPRA 0.3679233 0.006  0.04609683  0.06692226 POAPRA
## RANABO 0.3651906 0.007  0.06331065  0.06239837 RANABO
## RATPIN 0.7505251 0.001  0.29535902 -0.19033003 RATPIN
## SORNUT 0.7078727 0.001  0.42335123 -0.01950158 SORNUT
## TRIREP 0.4598473 0.002 -0.37948940 -0.03270281 TRIREP

2022 NMDS

adonis2(inner.hel.2022 ~ Treatment, data = env.data.2022,  dist="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = inner.hel.2022 ~ Treatment, data = env.data.2022, dist = "bray")
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  3  0.76026 0.25812 2.3195  0.001 ***
## Residual  20  2.18516 0.74188                  
## Total     23  2.94542 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Pair-wise comparisons
library(pairwiseAdonis)
## Loading required package: cluster
pairwise.adonis2(inner.hel.2022 ~ Treatment, data = env.data.2022, p.adjust = "bonferroni")
## $parent_call
## [1] "inner.hel.2022 ~ Treatment , strata = Null , permutations 999"
## 
## $LE_vs_LL
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.32134 0.24719 3.2835  0.004 **
## Residual  10  0.97863 0.75281                 
## Total     11  1.29997 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  1  0.08466 0.06313 0.6738   0.79
## Residual  10  1.25652 0.93687              
## Total     11  1.34118 1.00000              
## 
## $LE_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)  
## Treatment  1  0.31377 0.21079 2.6708  0.012 *
## Residual  10  1.17481 0.78921                
## Total     11  1.48859 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.24962 0.19812 2.4707  0.006 **
## Residual  10  1.01035 0.80188                 
## Total     11  1.25997 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)  
## Treatment  1  0.17717 0.16022 1.9078  0.037 *
## Residual  10  0.92864 0.83978                
## Total     11  1.10581 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.37395 0.23661 3.0994  0.009 **
## Residual  10  1.20653 0.76339                 
## Total     11  1.58048 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Year: 2023

## Call: capscale(formula = inner.hel.2023 ~ Treatment, data =
## env.data.2023, distance = "bray", add = TRUE)
## 
##               Inertia Proportion Rank
## Total          4.9591     1.0000     
## Constrained    1.3768     0.2776    3
## Unconstrained  3.5823     0.7224   20
## Inertia is Lingoes adjusted squared Bray distance 
## Species scores projected from 'inner.hel.2023' 
## 
## Eigenvalues for constrained axes:
##   CAP1   CAP2   CAP3 
## 0.7696 0.4783 0.1289 
## 
## Eigenvalues for unconstrained axes:
##   MDS1   MDS2   MDS3   MDS4   MDS5   MDS6   MDS7   MDS8 
## 0.7364 0.4885 0.3927 0.3212 0.2513 0.2101 0.1807 0.1551 
## (Showing 8 of 20 unconstrained eigenvalues)
## 
## Constant added to distances: 0.05150478
## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: capscale(formula = inner.hel.2023 ~ Treatment, data = env.data.2023, distance = "bray", add = TRUE)
##           Df SumOfSqs      F Pr(>F)    
## Treatment  3   1.3768 2.5622  0.001 ***
## Residual  20   3.5823                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Multiple comparisons for capscale for all contrasts of Treatment 
## 
## Model:  BiodiversityR::multiconstrained(method = "capscale", formula = inner.hel.2023 ~ Treatment, data = env.data.2023, distance = "bray", add = TRUE) 
## 
##             Df SumOfSqs      F Pr(>F)   
## LE vs. LL    1  0.56594 4.8862  0.002 **
## LE vs. NAT   1  0.22834 1.5823  0.098 . 
## LE vs. SIM   1  0.45445 3.3335  0.003 **
## LL vs. NAT   1  0.34988 2.7638  0.004 **
## LL vs. SIM   1  0.28030 2.5274  0.004 **
## NAT vs. SIM  1  0.59997 3.6140  0.003 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   UniqueBlock Block Treatment Year Bare_cover  Above Ground Middle      axis1
## 1        LE_1     1        LE 2023          3 1852.0  114.0 1305.5 -0.1492502
## 2        LL_1     1        LL 2023          3 1840.0   59.0 1363.0  0.5871939
## 3       NAT_1     1       NAT 2023          2 1854.5  127.0 1335.0 -0.6337099
## 4       SIM_1     1       SIM 2023         13 1819.5   87.5 1486.5  1.0284389
## 5        LE_2     2        LE 2023         18 1786.0  196.5 1252.0 -0.8450437
## 6        LL_2     2        LL 2023         46 1731.0  109.5  248.0  0.3190052
##        axis2 labels
## 1  0.8082363   LE_1
## 2 -1.1020888   LL_1
## 3 -0.6006288  NAT_1
## 4  0.9438200  SIM_1
## 5  0.8878608   LE_2
## 6 -0.4246796   LL_2
##                axis1        axis2 labels
## ACHMIL -0.0711105232  0.592520826 ACHMIL
## AGRHYE  0.0085664677  0.050758260 AGRHYE
## ALLVIN -0.0185623066 -0.016238688 ALLVIN
## AMBART -0.1280267347  0.032816374 AMBART
## AMBTRI -0.1574108999 -0.300615993 AMBTRI
## BARVUL -0.0269443322 -0.077481104 BARVUL
## BIDARI  0.0033968330  0.028002451 BIDARI
## BROJAP -0.0633261442 -0.103598384 BROJAP
## CARSPP -0.1212440666 -0.161028841 CARSPP
## CERSPP  0.0181368471  0.018000893 CERSPP
## CHAFAS  0.1606510797  0.055077603 CHAFAS
## CONCAN -0.0277846622 -0.052836330 CONCAN
## CORLAN -0.0212794604  0.241260701 CORLAN
## CORTRI  0.4820760468  0.112414607 CORTRI
## DESSPP -0.0157311613 -0.013651272 DESSPP
## ECHCRU  0.0223030327  0.019201888 ECHCRU
## ERASPE  0.0216923651 -0.019686356 ERASPE
## ERIANN -0.0132088157 -0.043483443 ERIANN
## ERYYUC  0.1372907701  0.002771190 ERYYUC
## FESARU -0.0135022460 -0.002617965 FESARU
## GALAPA -0.0146881661  0.013442529 GALAPA
## GALPED  0.0224825505  0.019236293 GALPED
## GERCAR  0.0108559078 -0.010491515 GERCAR
## GEUCAN -0.0142451724 -0.012361753 GEUCAN
## HELAUT -0.0583436417 -0.086514015 HELAUT
## HELMOL  0.3553255256 -0.034761103 HELMOL
## HORPUS -0.1475521122 -0.293480395 HORPUS
## IPOLAC -0.0090387556 -0.031331641 IPOLAC
## IVAANN  0.0455330518  0.039201868 IVAANN
## JUNINT  0.0136331062 -0.016777789 JUNINT
## KOEMAC -0.0487436827  0.170858654 KOEMAC
## LEPVIR -0.0572340552  0.003521346 LEPVIR
## LESCAP  0.1504633256 -0.021946339 LESCAP
## LIAPYC  0.0179812798  0.015481057 LIAPYC
## LONMAC  0.0004196733  0.031553343 LONMAC
## LOTCOR -0.0022111948  0.443190652 LOTCOR
## MONFIS -0.0331070975  0.064020687 MONFIS
## MYOVER -0.0378162987  0.046423163 MYOVER
## OXADIL -0.0721611624  0.020946220 OXADIL
## PANDIC -0.0207573701  0.008207711 PANDIC
## PENDIG -0.1192823538  0.025241468 PENDIG
## PHAARU -0.0859544746 -0.074776638 PHAARU
## PLALAN  0.0202676327  0.017207299 PLALAN
## PLAMAJ -0.0162510996  0.014872917 PLAMAJ
## PLAVIR -0.0141967833  0.012992818 PLAVIR
## POAPRA  0.0492069968  0.092727229 POAPRA
## PYCTEN -0.0189514955 -0.016445831 PYCTEN
## RANABO  0.0487820743  0.041907584 RANABO
## RATPIN  0.5240609418 -0.268029972 RATPIN
## RUDHIR -0.0093949298 -0.096825923 RUDHIR
## RUMCRI -0.0807045564  0.043057146 RUMCRI
## SETFAB -0.0327750423 -0.097865804 SETFAB
## SETPUM -0.4458985710  0.039650195 SETPUM
## SIDSPI -0.0161814285 -0.014042007 SIDSPI
## SOLCAR -0.0193228840 -0.027701254 SOLCAR
## SORNUT  0.7621760713  0.011031048 SORNUT
## SPHOBT  0.0372025018  0.032029647 SPHOBT
## SPOHET  0.0440449605  0.037920690 SPOHET
## SYMPIL -0.0346702003  0.027790858 SYMPIL
## TAROFF -0.0852401230 -0.004968236 TAROFF
## TRIPRA  0.0672797885  0.057924811 TRIPRA
## TRIREP -0.1045058986  0.087732099 TRIREP
## VERARV  0.0221452014  0.016543448 VERARV
## VERURT  0.0136331062 -0.016777789 VERURT
## VIOPED -0.0363283830  0.033247537 VIOPED
##   axis     ggplot        label
## 1    1 xlab.label CAP1 (15.5%)
## 2    2 ylab.label  CAP2 (9.6%)
##                  r     p         axis1        axis2 labels
## ACHMIL 0.557034573 0.001 -0.0711105232  0.592520826 ACHMIL
## AGRHYE 0.043162268 0.634  0.0085664677  0.050758260 AGRHYE
## ALLVIN 0.042659805 0.627 -0.0185623066 -0.016238688 ALLVIN
## AMBART 0.065102575 0.492 -0.1280267347  0.032816374 AMBART
## AMBTRI 0.332015883 0.018 -0.1574108999 -0.300615993 AMBTRI
## BARVUL 0.283979345 0.016 -0.0269443322 -0.077481104 BARVUL
## BIDARI 0.003289345 0.964  0.0033968330  0.028002451 BIDARI
## BROJAP 0.027646890 0.743 -0.0633261442 -0.103598384 BROJAP
## CARSPP 0.121420182 0.245 -0.1212440666 -0.161028841 CARSPP
## CERSPP 0.009249716 0.908  0.0181368471  0.018000893 CERSPP
## CHAFAS 0.287760954 0.019  0.1606510797  0.055077603 CHAFAS
## CONCAN 0.147106291 0.173 -0.0277846622 -0.052836330 CONCAN
## CORLAN 0.364442071 0.014 -0.0212794604  0.241260701 CORLAN
## CORTRI 0.635936293 0.001  0.4820760468  0.112414607 CORTRI
## DESSPP 0.027217602 0.918 -0.0157311613 -0.013651272 DESSPP
## ECHCRU 0.032565663 0.892  0.0223030327  0.019201888 ECHCRU
## ERASPE 0.045055339 0.611  0.0216923651 -0.019686356 ERASPE
## ERIANN 0.084905162 0.402 -0.0132088157 -0.043483443 ERIANN
## ERYYUC 0.211941307 0.092  0.1372907701  0.002771190 ERYYUC
## FESARU 0.017652522 0.812 -0.0135022460 -0.002617965 FESARU
## GALAPA 0.077753255 0.581 -0.0146881661  0.013442529 GALAPA
## GALPED 0.100065668 0.350  0.0224825505  0.019236293 GALPED
## GERCAR 0.009393616 0.898  0.0108559078 -0.010491515 GERCAR
## GEUCAN 0.022772173 1.000 -0.0142451724 -0.012361753 GEUCAN
## HELAUT 0.056859246 0.670 -0.0583436417 -0.086514015 HELAUT
## HELMOL 0.444901394 0.002  0.3553255256 -0.034761103 HELMOL
## HORPUS 0.312087172 0.017 -0.1475521122 -0.293480395 HORPUS
## IPOLAC 0.009490792 0.913 -0.0090387556 -0.031331641 IPOLAC
## IVAANN 0.311525948 0.008  0.0455330518  0.039201868 IVAANN
## JUNINT 0.133690558 0.202  0.0136331062 -0.016777789 JUNINT
## KOEMAC 0.370585711 0.009 -0.0487436827  0.170858654 KOEMAC
## LEPVIR 0.325601019 0.016 -0.0572340552  0.003521346 LEPVIR
## LESCAP 0.334806572 0.019  0.1504633256 -0.021946339 LESCAP
## LIAPYC 0.163347246 0.049  0.0179812798  0.015481057 LIAPYC
## LONMAC 0.172667641 0.105  0.0004196733  0.031553343 LONMAC
## LOTCOR 0.348615977 0.013 -0.0022111948  0.443190652 LOTCOR
## MONFIS 0.091609045 0.361 -0.0331070975  0.064020687 MONFIS
## MYOVER 0.173081368 0.121 -0.0378162987  0.046423163 MYOVER
## OXADIL 0.143002494 0.209 -0.0721611624  0.020946220 OXADIL
## PANDIC 0.006977944 0.921 -0.0207573701  0.008207711 PANDIC
## PENDIG 0.193114410 0.108 -0.1192823538  0.025241468 PENDIG
## PHAARU 0.041539879 0.794 -0.0859544746 -0.074776638 PHAARU
## PLALAN 0.001264554 0.986  0.0202676327  0.017207299 PLALAN
## PLAMAJ 0.132069256 0.231 -0.0162510996  0.014872917 PLAMAJ
## PLAVIR 0.091934859 0.479 -0.0141967833  0.012992818 PLAVIR
## POAPRA 0.137217443 0.219  0.0492069968  0.092727229 POAPRA
## PYCTEN 0.083336781 0.521 -0.0189514955 -0.016445831 PYCTEN
## RANABO 0.181282020 0.125  0.0487820743  0.041907584 RANABO
## RATPIN 0.671445336 0.001  0.5240609418 -0.268029972 RATPIN
## RUDHIR 0.049085637 0.584 -0.0093949298 -0.096825923 RUDHIR
## RUMCRI 0.068456245 0.460 -0.0807045564  0.043057146 RUMCRI
## SETFAB 0.118809057 0.255 -0.0327750423 -0.097865804 SETFAB
## SETPUM 0.545215500 0.001 -0.4458985710  0.039650195 SETPUM
## SIDSPI 0.059385684 0.677 -0.0161814285 -0.014042007 SIDSPI
## SOLCAR 0.007758942 0.916 -0.0193228840 -0.027701254 SOLCAR
## SORNUT 0.789440722 0.001  0.7621760713  0.011031048 SORNUT
## SPHOBT 0.171127596 0.118  0.0372025018  0.032029647 SPHOBT
## SPOHET 0.163347246 0.049  0.0440449605  0.037920690 SPOHET
## SYMPIL 0.105918848 0.330 -0.0346702003  0.027790858 SYMPIL
## TAROFF 0.111769896 0.291 -0.0852401230 -0.004968236 TAROFF
## TRIPRA 0.163347246 0.049  0.0672797885  0.057924811 TRIPRA
## TRIREP 0.412564825 0.003 -0.1045058986  0.087732099 TRIREP
## VERARV 0.005195253 0.958  0.0221452014  0.016543448 VERARV
## VERURT 0.133690558 0.202  0.0136331062 -0.016777789 VERURT
## VIOPED 0.211880088 0.036 -0.0363283830  0.033247537 VIOPED
##                r     p       axis1       axis2 labels
## ACHMIL 0.5570346 0.001 -0.07111052  0.59252083 ACHMIL
## CORTRI 0.6359363 0.001  0.48207605  0.11241461 CORTRI
## HELMOL 0.4449014 0.002  0.35532553 -0.03476110 HELMOL
## IVAANN 0.3115259 0.008  0.04553305  0.03920187 IVAANN
## KOEMAC 0.3705857 0.009 -0.04874368  0.17085865 KOEMAC
## RATPIN 0.6714453 0.001  0.52406094 -0.26802997 RATPIN
## SETPUM 0.5452155 0.001 -0.44589857  0.03965020 SETPUM
## SORNUT 0.7894407 0.001  0.76217607  0.01103105 SORNUT
## TRIREP 0.4125648 0.003 -0.10450590  0.08773210 TRIREP

2023 NMDS

adonis2(inner.hel.2023 ~ Treatment, data = env.data.2023,  dist="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = inner.hel.2023 ~ Treatment, data = env.data.2023, dist = "bray")
##           Df SumOfSqs      R2      F Pr(>F)    
## Treatment  3   1.2223 0.32383 3.1928  0.001 ***
## Residual  20   2.5522 0.67617                  
## Total     23   3.7745 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Pair-wise comparisons
library(pairwiseAdonis)
pairwise.adonis2(inner.hel.2023 ~ Treatment, data = env.data.2023, p.adjust = "bonferroni")
## $parent_call
## [1] "inner.hel.2023 ~ Treatment , strata = Null , permutations 999"
## 
## $LE_vs_LL
##           Df SumOfSqs     R2      F Pr(>F)   
## Treatment  1  0.55758 0.3416 5.1883  0.004 **
## Residual  10  1.07469 0.6584                 
## Total     11  1.63228 1.0000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LE_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)
## Treatment  1  0.22834 0.13662 1.5823  0.106
## Residual  10  1.44310 0.86338              
## Total     11  1.67144 1.00000              
## 
## $LE_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.44745 0.25705 3.4599  0.003 **
## Residual  10  1.29325 0.74295                 
## Total     11  1.74070 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_NAT
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.34918 0.21714 2.7736  0.005 **
## Residual  10  1.25892 0.78286                 
## Total     11  1.60810 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $LL_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1   0.2803 0.20175 2.5274  0.004 **
## Residual  10   1.1091 0.79825                 
## Total     11   1.3894 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $NAT_vs_SIM
##           Df SumOfSqs      R2      F Pr(>F)   
## Treatment  1  0.58171 0.28249 3.9372  0.006 **
## Residual  10  1.47748 0.71751                 
## Total     11  2.05919 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Plot Panel

Question 3b: Does timing of species arrival influence species composition? - Jaccard

Year: 2022

Year: 2023

Plot Panel